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I. INTRODUCTION 


This work is the continuation and the end of a cooperative effort mad? by the following 
persons: 

— ► Hieu HAMINH, Professor at INP, Toulouse, France, 

— ► Dany VANDROMME, Professor at INS A, Rouen, France, 

— ► Wolfgang KOLLMANN, Professor at UC Davis, 

— * John VIEGAS, Research scientist at ARC, 

— + Morris RUBESIN, Senior Research scientist at ARC. 

to develop a second order closure turbulence model for compressible flows and to implement 
it in a 2D Reynolds- averaged Navier-Stokes solver. This work has been initiated in the 
early 80’s while one of the authors (DV) was NRC Research Associate with the Ames CFD 
branch. During the subsequent years, the NASA grant NCC2-186 allowed the continuation 
of this work through repeated funding from the Experimental Fluid Dynamics branch, lead 
by Dr. Joseph MARVIN. 

From the beginning of this work where a k — t turbulence model was implemented 
in the bidiagonal implicit method of MACCORMACK (referred to as the MAC3 code) to 
the final stage of implementing a full second order closure in the efficient line Gauss-Seidel 
algorithm, numerous work have been done, individually and collectively by the individuals 
mentionned above. “ 

Besides the collaboration itself, the final product of this work is a second order closure 
derived from the Launder, Reece and Rodi model to account for near wall effects, which 
has been called FRAME model, which stands for FRench-AMerican-Effort. 

Another benefit of this collaboration was the proposition and extensive testing of 
various turbulence model corrections to account for strong compressibility effects. Among 
the various contributions in this field, the following main lines has been worked out: 

-+ The modelling of the pressure and density correlations based on, among other as- 
sumptions, the polytropic assumption. This approach has been initiated early in the 
70’s by Rubesin, and taken over by Vandromme. 

— ► Kollmann and Vandromme have introduced the compressible version of the e equation 
with specific compressibility corrections mostly based on the mean velocity divergence. 

— ► Later, the various proposals based on the compressible dissipation made independently 
by Sarkar and Zeman has been tested also by Viegas and Rubesin and compared to 
the various Rubesin proposals for the compressible mixing layer. 

— ► More recently, Vandromme continued to work on new models for the pressure dilata- 
tion in presence of strong shocks. This work, which has been conducted during a work 
at the Center for Turbulence Research with Zeman aimed also to cross-check earlier 
assumptions by Rubesin and Vandromme. 

In common with all the contributions which have been done under the ARC-UCD 
grant, the authors must recognize that they spent a lot of time to play with the numer- 
ics, rather with the strict model of physics of turbulence. That confirmed, as a general 
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conclusion that: 

-+ a turbulence model is never independent of the numerics, 

— + there is no hope to develop a universal model of turbulence, the best to be expected 
being a well-calibrated model for a given type of flow, 

— > a turbulence can not be used as a black box without a minimum of expertise 


During this last summer period, two different problems have been worked out. The 
first was to provide Ames researchers with a reliable compressible boundary layer code 
including a wide collection of turbulence models for quick testing of new terms, both in 
two equations (Jfc - e, k - w 2 , q - w, k - kl etc. . . ) and in second order closure (LRR and 
FRAME). The second topic was to complete the implementation of the FRAME model 
in the MAC5 code. The work related to these two different contributions is reported in 
the following two chapters. 


II. NUMERICAL PROCEDURE FOR BOUNDARY-LAYER EQUATIONS 


II-l. INTRODUCTION. 

Originated by PATANKAR and SPALDING [1], the numerical procedure presented here 
contains several specific modifications, concerning more particularly 
i) the treatment of source terms, 
it) the slip false grid points, 

Hi) the boundary conditions, and 
iv) the numerical algorithm 

used to solve the coupled partied derivative parabolic equations. 

The procedure concerns a set of partial derivative transport equations including: 

• Continuity equations: 

+ = o (II -» 


• n transport equations: 


U llc*" +V Tr*" 



(II - 2) 


II-2. TRANSFORMATION OF THE EQUATION SET 

To avoid the continuity equation (II-l), we can define a stream function xp as: 



(II - 3a) 
(II - 36) 
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Therefore, all variable function $(X, r) will be considered as function of the VON 
MISES variables x and ip — ► x being the longitudinal coordinate on the interned 

boundary of the flow field (Fig.l). 



* Remark: Neither the internal boundary I nor the external boundary E is necessary a 
streamiline, except if the equations are no longer parabolic. 

The elementary flow rate through the crown area ds = 2irr.dr is: 

dQ = 2tt. r. dr. p.U cos (d = 2Trcos(d. dip 

(/d(x) being the local angle of the velocity vector). 

Assuming that (d <C 1 — ► cos/? ~ 1: 

X = X(x,ip) r = r(x,ip) (II — 4) 


with 



x = X 

dip = p.U.r.dr 


(—) =i 

\dXJr 

( dx \ 


\drJx 

/dip\ 
V dx > 

| = -pV r Q 

r 

(dip\ 

\dr)x 


i i d r „ n (d$ 


J. 
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Therefore, the equation (H-2) becomes: 



(II - 5) 
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Let us introduce now new variables defined as: 


x => x = x 


I ^ _ t ~ 


u is function of ip and x (through ipi and iPe) 

(=)-» 


(— ) = 

\dib) 


dip / ipE — ipi 


(du\ 

ip - 

- ip I 

dipE 

dip r 

1 dipi 

vae/. 

(ipE - 

~ *P I ) 2 

di 

d£ J 

ipE — ipi 


a$ dx * du 

\ dx ) ip dx* dx ' du dx 


with: 

Defining: 


/ 5$ > 

1 

[M +W l 

f dip e 

dip j\1 

wx* > 

'w (ipE-lpl) 

la* +w| 

\ dx 

/J 


a$ dx* d$ du 
\dipy x dx* dip du dip 

1 /d$\ 

( IpE ~ ipi) ' du ) x* 


a 

r r a$i 

1 a r 

K a$i 

a$i 

5t/> 


(ipE — Ipi) du L 

1 a r 

^av*. 

AC 


1 ^ 

II 1 

(ipE ~ 

■ a$ 

ipi ) aw J 

(V>E - V ’/) 2 

. aw 


K = r 2a p.U.V 


rhl = pVi = . 


1 dipi 


dx 


T/ 1 
m* Vs = 


a$ > 

i i 1 - f 

lax*> 

(V’E - V 5 /) ■ 


'd$\ 

<du ) i 


r , a$l S* 


(ipE ~ V’/) 2 




/ 5$ \ . , , /a$\ 5 / a$\ 

(dd„ +(a+M (d,- = d c d +d 


a = 


r/m/ 
1p E - ip I 


b = 


r_Ern£ - rjrhj_ 
IpE ~ ip I 


C = 


pUf^D 

(rpE ~ ipl) 2 


d = 


S* 

~pU 


(II - 6) 


(// - 7) 


(77 - 8) 


(77 - 9) 


(77 - 10) 


(77 - 12) 
(77 - 13) 
(77 - 14) 
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II-3. ORIGINAL EQUATIONS TO BE SOLVED 


± { ,u )+ f y m=o 


u^ + vf 

dx dy 
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P uv dU 
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( n - 1 )C P ' T 2 ' dy 
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dx 
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dy 
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u<hal + v— 

dx^ dy 
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(// - 21) 


(// - 22) 


In the successive transformations: 

(x,y) = 
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Therefore, the general form for all equations to be solved is now written as: 


r rrd$k 


U 


(7/ - 23) 


dx 


A a f a$n 

M -- 

/=1 

M 

+ ^ Rki$i + Qk 


A = 

Cki = 


mj 


ipE - ipl 

pV£ki 

(rp E - ip!) 2 

At/ 


Ru = 


P-U 


1=1 


B = 


THE — TUI 

&kl 


Dki = 


XpE ~ V’i 


Qk ~ p.U 


(II - 24) 


II-4.- THE BLOCK-SOLVER TECHNIQUE. 

The previous P.D.E. could be solved by using a Block-Solver technique (See Ref.[2]). 
This technique is presented as following: 

Denote <J>; = and integrate over V; (Fig.2): 


Figure 2 
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d$ 

Variation of -r— with u > : 

ox 



/d£\ _ ~ $7 

\ dx / i Ax Ax 


Hence: 


^~T)SL dxdu - f = (f)- 


Ax(wj + i 


= -( / / do;*") 

Ai(wi+i -Wj_i) J > 

-l [ ' du$+ f ' + ' du$- f du;$--[ 
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i + l 
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Finally: 
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Wi - w tl 


+ 3 + $i+i 
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" — - at- - *r+. 1 " i+l ” 
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<^+l — Q7j 1 
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(II - 25) 
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d$ 

• Discretization of: with A = constant : 

Ou> 


2 / f dx.du.A.— = (.4.—) 

Ai(w i+1 -ui-i)JJ Vi du \ du>Ji 

/ A d$\ 1 . r p , A d$ P+ l , 

( A ^)rvML^ + L 


cL> 1 


/ , d$\ 2 A f f Wi , p +l , 1 
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Finally: 


U.£) . A . 

V do; / > u>i+i — w,_i 


(/J - 26) 


• Discretization of: B.u> 


du' 


with B = constant: 


1 

V, 



dx.dw.B. 


dPhi 

du> 




2 B 


^i+l ~W»-1 



5$ , 

uj—du + 

atjj 



Finally: 


r> <&\ B ( * 3w, +wi_i , ^ 3u>i + Wi+i } fTT ^ 

(B.W.— ) = + > (11-27) 

\ ou / i 4 1 u;,-+i — u; : _i a; l+ i — u;^ J 


/ (7<p \ 

Discretization of: (C. —J : 
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« (w l+1 
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Finally: 
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[(<?£) 1 =7 ‘ r{ + *,-,(^^) 

L\ du/ii (Wi+l — W|_l) l \Ui — LJi-\ / 


(wi+ 1 


^ (Ci + Ci+i ' Ci + Ci- 1 

1 

'^i + 1 “Wj 

(O + CW,)} 

Vwj+i — Ui / i 


+ 


(// - 28 ) 


5 $ 

• Discretization of: Di(u>)-t—: 

du> 


- (»S), 




(II - 29) 


Finally: 


( D l), - + 31>i)( *' - *'-•) + (I)i+ ' + 3A)(4i+1 - *■>} 


The finite-difference equation. 
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Ui ~ Ui - i - + 3 *‘ + »? +1 — 
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i-^i 1 
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W,+l — Wj_i 
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1 

+- 
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M 
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M 


+ '£ R ‘ i ‘ + Q < 


1=1 


(II - 30) 


k = 1,2, ...M 


By multiplying by (u>i+i — Wi-i), we obtain: 


4 Ax l ^* 7 * — + 3 ( w^ + i — + ( w<+i — w *)$ f +1 

— (wj — + 3(u>j+i — + (wj+i — 
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Otherwise: 


, 1 B T'-'i-l 

+ ^f-n — 77 A ~ t( 3u, « +^«-i) ' 

*1 4Ax 4 u>i—Ui - 1 


C kk + Cgi + 3D* 


Jtfc 

-} 


+**( 
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4Ax 
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.t x , fl,, , . C,“ + C“ D&+3 Kf‘ 

+ **+■( 4 A ~ +/ 1+ 7 ( ^ i+ " | -‘ ) ~ 


u>j+l — Wi 


+ 


M 

£*!->{ 


0‘ii+3Cf' cf* + c“. 




or + o,_i i 

U>j — CJ,_1 J 


4- 


f *- { c - * + c “' + c " - + c “-. + 


M 


E *!♦.{- 


Dfi i + 37)*' C** + C kk 


Wk 

= Qk(ui + 1 — u;,_i) + 


4 w,-+i 

w,- -w«-i ,_jt 


u »+i 1 
— U?i J 


4A* *- + JE'"** - W - )<tr ‘ + 


4Ax "‘ +1 

(77 - 32) 


Let us call: 


Zb = 


— A — — (w,+i 4- 3o? t ) 
B 

+ A + — + 3a >i) 
4Ax 4 

B' 


UJi+l -Uj 

4Ax 

- w,_i 


(77 - 33) 


Zc = (lAl ■ 7) (u,i -' ~“ i+,) 
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To apply a block-solver algorithm, let us put this system into this form: 

+ A i.xi + Bi.Xi+i = yi (II - 34) 

with 2 < i < N\ . Except: 

Ai.xi + B 1 .X 2 + C 1 .X 3 = yi (II — 35) 

Bn-xn-2 + Ch.xn 1 + An-xn = VN (11-36) 

In ( 34 ), ( 35 ), (36) A,,Bi,C, bring (M x M) matrices, Xi and yi being (M) vectors 
(M=number of equations to be solved). 

Then, the matrices A*,Bi and C, are: 


A = 


_ c? + c*U + cj“ + eft 1 _ DfU - ) 

. cjj-j-i — cj| u?i — uj 1 — 1 4 

+ - Wj_] ) - 2a for l = k 

c? + c‘l, + C* 1 + eg, 1 _ a” 1 - of'. 

L Wi+l — U>i U>i —OJi-1 J 4 

{ + R^Ui+i -U7j_ 1 ) for l ^ k 


(II - 37) 


B,* = 


, Cf + CtU , A>i + 37>f' , „ 

-| 1 1- £j A 

Wj+1 -u>i 4 

+ cf + c !L + hl/l 

wj+i — Wj 4 


for / = it 


(77 - 38) 


cJ 


C ? 1 + C& 7?fl, + 37)*' 


+ — — + Z B for / = fc 

+ for,/, 

u;,- — u>i_i 4 


yf = — 


1 


(Wj — W,_i)$j_ fc i + 3(u?,+ l — 0!i_i)$j * + (Wj+i ~LOi)$i + 1 


-k 




4Ax 

- (wf+i -ui-\).Qk 


(II - 39) 
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II-5.- BOUNDARY CONDITIONS. 


The previous relationships are only valid for all points i from 3 to (N - 2). We have to 
modify the previous results for the point 2, assuming an integration to the wall, and for 
the point (N - 1), assuming an integration to the outer boundary. 


A* = 


C\ x + Cf* 3 C? + Cj‘ 


Ll>2 — CJl 2 (u?3 — U)2 J 

4 £>*' + £>*' - 

*: h R 2 ( w 3 + ^2 — 2o;x ) 


-A — f - SaJgLifa i 


(II - 40) 


B k 2 = 


3 C kl 4- C kl B v w 3—0*2 

A — —(0*3 + 3cj2 ) — 

u>2 — u\ 4 4Ax 

3 D\ l + D kl 


(II - 41) 




C ‘“ + C " +2 A- B( Ul + Ut ) - + flf')) 


U)2 ~ U>1 


(II - 42) 


y 2 =- 


(0*3 — U)2 )^3 * + (3tL>3 + U>2 — 4u*i )$ 2 * + 4(u>2 — U\ )$ j 


k 


4Ax 

— (u$ + 0*2 ~ 2u*l )-^2 


(II - 43) 


A^ — < 


Cjy + Cjy-i , 3C^ + C^/ 2 


+ 


ojn—un- 1 2(u;Ar t - WjVi J 

40ft + 0ft, - 0ft 


B 

+ A + — (4o>at + w/v, — u>;v 2 )- 


, D *J /o 4a; /v - W/Vj - 3o*/v 2 

+ UnA^N -uni -un,) — - 

4Zxr 


(II - 44) 


C k - 
^N x ~ 


3 C Nx+ C N, 1 \ UNx-UN, 

+ iST" 

3 pfj + *>y, 

4 


(II - 45) 


Bft = | 


eft + g& 

Wat — 0> J V_ 1 


2A — B(un + ) — 


un — un 


Ax 


HK + 0ft) } 


(II - 46) 
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y k Ni = 


4Ax 


[(WN-I - UN-2)$N-2 + -WJV, - <$WN 7 )$N k -i +4 (un -<*>N- l)$/v* 


- (2u>N —UNi - UN 2 )-Qn-2 


(II - 47) 
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III.- IMPLEMENTATION OF THE COMPRESSIBLE SECOND ORDER 
CLOSURE IN THE MAC5 CODE 

From the various attempts which have been worked out previously, a simple method is to 
be used for the treatment of the source term collection, to ease their numerical treatment. 


III-l.- The FRAME model 

The second order moment closures (or Reynolds stress models) are currently the most 
general one-point correlation models from the point of view of physical theory. These 
models require the solution of additional field equations for the complete set of Reynolds 
stresses the turbulent heat flux vector pv a T” and, frequently a scale equation 

which can be eorw 2 , similar to that one used for the two-equation models. These models 
axe obviously more complicated than the eddy viscosity based models. One of the most 
important physical properties contained in these models is however a stress relaxation 
property which cannot be correctly represented in the eddy viscosity models. 

From a simple manipulation of the instantaneous Navier-Stokes equations, derivation 
of a transport equation for the Reynolds stress components yields: 


j-(p ) + j£-\p v m a vjf V. y + P vlvpv* + s a 7 v 0 p' + v Q p' - p{S Q1 Vp + S M v a ) 

dvl dv] 


_ r— dv 0 dv a -\ ,.(dv Q ovp\ 

= -f aZ + ” ^ azj + p \dZ + dTJ 


P 


dVfi dvli 

7 


/Q7 ax. ' 


— dp — dp 
V °dxp V ^dx Q 

{III - 1) 

By taking advantage of the contracted index convention (a = /? and summation) and 
dividing the resulting equation by 2, we obtain the exact form of the transport equation 
for the turbulent kinetic energy: 


d . 


(p l~, k + pv"k + v'' y p’ - p »’5„,) 


_ n n dv a 

- p v ^dZ 


+ 


dx 


dx a 


,dv a -irdp dv a 


(III - 2) 


The modelling of the turbulent stresses and fluxes introduces a lot of new terms, so far 
not defined by experiments. It is not possible at this time to reach definitive conclusions 
on the validity of all these closures. For this reason, in the following, emphasis is put 
on the modelling of the Reynolds stresses whereas the remaining turbulent fluxes will be 
handled with a general anisotropic form of the gradient approximation. 

This remark seems very restrictive, but in fact the lack of experimental results makes 
the Reynolds stress modelling problem sufficiently complex to delay the equivalent treat- 
ment of others fluxes. Practically, the application domain will be restricted to compressible 
flows with moderate heat and mass transfers, although extension for combustion flows has 
been made already [3]. 
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Let us return now to the open transport equation for pv" a v"p (eq. (Ill- 1)) in which all 
new unknown terms need to be modelled. 

This equation can be modelled by extending the incompressible models of Launder, 
Reece and Rodi [4] and Hanjalic and Launder [5] to compressible conditions, i.e. using 
Favre decomposition, introducing the non-zero divergence terms that were eliminated in the 
original models and accounting for the non- zero mean mass-weighted fluctuating velocities 
[6]. In this report, we restrict the discussion to the important points of modelling. 

• turbulent flux of Reynolds stress : = pu” 

Starting from the exact transport equation for pv” a v”gV y , it is possible, as shown by 
Hanjalic and Launder [5], by neglecting diffusive and convective terms, to obtain the 
following form: 


— ” ” ” ^ ^ » » n n d nil, » » ^ »»\ / r V r A \ 

-p v Q VpV y = C a p ~(v Q v 6 —v p v y + v 0 v 6 —v a v y + Vy v 6 —v Q v 0 ) (III- 3) 


This form conserves the symmetric character of the third order tensor but, for 

practical purposes, a simpler form, suggested by Daly and Harlow [7] seems to produce 
results of similar quality. 


_ n » n />( _ k „ » d n n 

-w-Vi = C 'P- v i v ‘e7, v « v f 


(III - 4) 


• pressure diffusion: = v" a p' 6p y + v" 0 p' S ay 

Most people neglect pressure-induced diffusion term, mainly due to the lack of exper- 
imental information. The measurements of Irwin [8] in a wall jet suggest that this term 
cannot be very important. Furthermore, some authors argued that the pressure induced 
diffusion, if non negligible, would act to destroy the symmetry character of the triple 
correlation term and support the use of the compact form given by equation (III-4). 


• viscous diffusion term: = 


dx. 


( v a P Spt v p P 

Assuming that a) the correlation between viscosity fluctuations and other quantities 
is weak, b) the product of density correlations with velocity gradients is small. Then the 
development of the molecular diffusion term is written as: 


. d ( d , d / „ dv'p »dv"\ 

v,sc. <UF. = —^— Va ^+—^ Va — +tlv ^ 


If the flow is incompressible or solenoidal (weak compressibility, Dussauge [9]), the 
viscous diffusion can be written as: 


visc . diff . = a ( ^ 


(III - 5) 
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* pressure strain correlation: 


=?& + £i) 

' dxp dx a 


In strongly anisotropic flows, i.e. in situations where second order closures are needed, 
this term is a central piece for explaining the redistribution mechanism between Reynolds 
stresses. To model this term, the approach is an incompressible-like technique, which 
consists in integrating a Poisson equation for the fluctuating pressure. The result of this 
integration, transposed to variable density flows is written as: 




C*2 4- 8 
11 





where 


and 


~ 30< 55 ~P kS -> ~ - f 4 ««*) 

+ 7rS c ^^ “ I 4 "'’** + C,( ^ - 

_ , ,dv^ dva 2 C dv-y,\ 

+C5M( a7^ + a^; _ 3 i “'W) 

(III - 6) 

p = D = -pvy 9 ^ = p„ 

(III - 7) 

CP’'" - / n dv a x 

P *> = -P K^+v^) 

(III - 8) 


(III - 9) 

, -2.5 

/l=eXp l + iW50 

(III - 10) 


The first two lines of equation (III-6) represent the redistribution mechanism in the 
flow field far from the wall (Launder, Reece and Rodi [4]) and the last two lines of this 
equation take into account the wall influence in this mechanism (Hanjalic and Launder 
[5]). The effect of this last contribution is twofold: 

* it has an opposite effect to the classical return to isotropy term of Rotta [10]. 

* it acts also as a rapid term to increase the anisotropy of the stress production terms. 

It must be emphasized that the transposition of an original incompressible technique 

to a variable density situation is not free of uncertainties. For instance a corrective term 
appears in the development of the non linear contribution due to the use of Favre averaging 
(see [11]). Also, the fourth rank tensor, corresponding to the high Reynolds number rapid 
term does not possess all the mathematical properties of its incompressible counterpart, 

b li ± b s2 or b s2 ± 0 

Finally, the whole term can be considered as a pure redistribution contribution only 
in the case of solenoidal turbulence field. Otherwise the bulk deformation is a source/sink 
term (Dussauge [9], Vandromme [12]). 
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The real weakness of this approach is the use of an incompressible approach (Poisson 
equation) when the flow is compressible. It would be more justified to introduce a wave- 
like operator to evaluate the fluctuating pressure field, like Feiereisen [13] tried from the 
results of his direct simulations. Unfortunately, the results, which have been obtained are 
nearly identical to those of the incompressible formulation. 


* mean pressure gradient term: = — v' a 

oxp 


— dp 
~ Vf, dx a 


According to Rubesin’s proposal for the one-equation turbulence model, that term 
can be treated like compressibility terms for two-equation turbulence models, assuming a 
polytropic behaviour of the fluid (see chapter 6). 


* viscous dissipation 

It has been shown, in incompressible case [14], that the dissipation tensor is diago- 
nally dominant and nearly spherical. The ratio of the deviatoric to the diagonal terms 
being related to the Reynolds stress anisotropy, the dissipation is described with a com- 
pound function which is scalar in the high turbulent Reynolds number zones and allows 
an anisotropic dissipation elsewhere (wall vicinity for instance). 



v a v 0 


/. + 0 - /.)§*«»*) 


(Ill -11) 


with 

fa = 1/(1 + Re t /10) ; Re t = k 2 /ue 

Nevertheless, some of the ’’slow” pressure strain terms may also represent anisotropic tij. 

To summarize the assumptions made above, the modelled Reynolds stress equation 
can be written as: 


d , _ \ , d / _ » n j-tt —k „ » d » n d 

j f (P v Q Vi J )+—{pv-rV a v l} - C a p-v y v 6 —v a v 0 - p—v a v fi ) 


dxs 


dx. 


r7 — ” C*2 + 8, -jr — ' 2 c pr' . 8C 2 2 . JT — 2 . ~ . 

30C2 — 2 _ n _ € ( T - ” 2 

^ P K S a p - Cip-{v a v 0 - -d a p k) 

+ v£: ) - + (1 ' 

■&} - | s a ,k) + c t (Ki- 




( III - 12) 


The last unknown, which remains in the modelled Reynolds stress transport equation, 
is the turbulent dissipation rate e. The modelling of a transport equation for this quantity 
has been given already for the first two equation model, and only the discrepancies due to 
the different level of closure are of some interest here. 


21 


A basic difference compared to the eddy viscosity model is that the Reynolds stresses 
can be considered now as exact quantities. This yields a more accurate evaluation of P*, 
the production term. Furthermore, as the eddy viscosity does not exist any longer, the 
turbulent diffusion transport of the dissipation is modelled by a generalized gradient flux 
approximation from Launder [15]: 


n k n ” 

-„ ae = C,- c V < ,v ll — 


(III - 13) 


As far as compressibility terms are concerned for the dissipation equation, the exact 
derivation of the equation introduces variable density terms as shown in [11] but the usual 
method is to ignore these terms and keep the e-equation similar to its incompressible 
counterpart. The modelling is done globally. If compressibility terms are introduced in 
the turbulent energy equation or here, in the Reynolds stress equations, experience shows 
that their counterpart is needed in the dissipation equation as well. 

The modifications induced for the total energy equation are derived similarly; the 
turbulent fluxes are expressed with an anisotropic relationship and the triple correlations 
follow the same approximation as in the Reynolds stress equation. The modelled total 
energy equation writes now as: 




/ / 2L 2 \ / _ >T n dvy . 

{{pv-r E + p + + KP v a v i - P-FTT) Ua 


dx. 
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ln n _k~. , dT , dv 7 dvp 2 dv a .-ir\ _ 

-(C. C,-l/,-V a V y +S fy h)j- - + g- - 5 6 ltdZ> U i>) - 


dv^ 




dx p 


dxp 5x 7 3 1 dx a 


treatment. 


III-2.- Implementation 

For many years, the development of numerical methods has been motivated by the need 
of solutions for the Navier-Stokes equations. Only recently, has interest increased in the 
solution of turbulence models and the development of accurate turbulence models has been 
recognized as a necessary route. Indeed, the interest in algebraic models has been due in 
part to their inherent simplicity, but also to the straightforward extension from laminar to 
turbulent cases by merely an alternate definition of the viscosity coefficient. Unfortunately, 
experience has shown that such a crude modelling assumption was not satisfactory as 
soon as the flow was slightly complex. The use of transport equation turbulence models 
introduces the turbulent kinetic energy, which needs to be accounted for in the total energy 
budget. For incompressible flows, this concept is not relevant, since the pressure is not a 
thermodynamic variable, but has only a mechanical role. In compressible flow the situation 
is quite different and the existence of k is felt everywhere in a Navier-Stokes solver, even 
inside the Euler part. 

A second difficulty, which is associated with transport equations for turbulence models 
is the treatment of non conservative source terms. As most of numerical schemes for Navier- 
Stokes equations took advantage of their strong conservative character, problems related 
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to the stability and the stiffness of source terms has often been discarded. We will examine 
in this report, various techniques to handle these problems, especially in the framework of 
implicit schemes. 


II-2-1. Energy coupling 


The instantaneous form of the total energy definition is: 


p E = p e + i pv Q v a (III - 14) 

In terms of Favre mean and fluctuating components, this equation becomes, after time 
averaging: 

E = e + ^v a v a + k (III -15) 

Therefore the solution of the temperature field from the total energy equation requires 
the knowledge of the turbulent kinetic energy. Neglecting that quantity [16] is equivalent 
to ignoring the energy which is extracted from the mean motion to constitute the turbu- 
lence energy. For incompressible flows, this is ignored and the turbulent motion is only 
superimposed to the mean. The coupling appears only through the mechanical role of the 
turbulent stresses, which are added to the viscous terms. For compressible flow calcula- 
tions, all the energy exchanges between the various scale motions must be considered to 
satisfy the global energy budget. It is well understood that, in most of the inviscid part of 
a flow field, the turbulence level is very low and the energy budget is not affected. But in 
regions with high shear or strong pressure gradients, the turbulent kinetic energy can be 
of the order of the mean, and must not be neglected as done usually [17], [18], [19]. 

The complete formulation of the constitutive relationship for the Reynolds stress is 
written in terms of density weighted variables as: 


_ - 7 — ». rdv a dv 0 2 dv y , 2 

-» v ° v > = te; + di: - ; ■ e °> srl - 5 ■ s °> r‘ 1 


Sly 


(III- 16) 


in which the turbulent kinetic energy term makes the contracted index form possible. This 
feature appears explicitly in the momentum and total energy equations where a turbulent 
normal stress is added to the mean pressure. A so-called effective pressure can be defined 
in the following way: 

p*=p + ^pk (III -17) 

In fact, this turbulent contribution to the pressure field is only an approximation 
neglecting the anisotropic nature of the Reynolds stress tensor. It was only introduced to 
insure a non-zero trace of this tensor. 

Such an approximation does not take place within the framework of a second order 
closure. In that case, the normal stresses appear explicitely in the momentum and total 
energy equations. Then, the relevant effective pressure is not isotropic any longer, but 
also depends on the turbulent energy distribution on its three normal components. For 
instance, in the o-momentum equation, the effective pressure will be: 


P* a =P + pv" a 2 


(III - 18) 
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Unfortunately, the concept of an anisotropic pressure field is difficult to handle, es- 
pecially with respect to the temperature field. Therefore, it is necessary to follow the 
same reasonning as for the static pressure definition from the kinetic theory of gases, and 
approximate the turbulent pressure by the mean of three components, i.e. | pk. 


III-2-2. Diagonalization of jacobian matrices 

To avoid the severe limitations of explicit methods, implicit schemes are preferred. A 
classical (but non unique) method for obtaining an implicit approximation is to take the 
time derivative of the original system. 


d_ 

dt 


dU dF dG = 
. dt dx dy 
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with the vector elements: 


U = 


' p ' 


t 

'll 
1 

pu 


pu 2 + p -\-JPF 2 

pv 


puv + pu" v" 

pJL 


puE + (p + pu” 2 )u + pu” v” V 

pu” v" 

F = 

puu”v” 

Pi 


put 

pxF 2 


put/” 2 

pv” 2 


puv” 2 

. pw ” 2 . 


puuF 2 


G = 


pv __ 
puv + puv" 
pv 2 + p_j- pv” 2 
pvE + (p + pv" 2 )v -|- pu” v" u 

pvu”v” 

pvt 

pvu ” 2 
put;” 2 
pv w” 2 


Define the jacobian matrices as: 


A - — ■ 

du ' 


B= dG 

dU 


r-^L 

6 dU 


(III - 20) 


the implicit approximation writes as: 


/ r . dA» dB • 

(/ + A '"aT + Ai 'a7 


A t.C)6U n+l = A U r 


(III - 21) 
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with the following increments: 


6U n+1 


. dU n+l 


AU n+1 



( III - 22) 


Equation (III-21) can be solved either by approximate factorization or by classical 
relaxation methods such as line Gauss-Seidel or point Jacobi. 

Equation (III-15) is used in the development of the diagonal form of the jacobian 
matrices 4 and B. Consider, for instance the x-direction, the jacobian A can be related 
to its diagonal form by the relation: 


A = 


( III - 23) 


To illustrate that, only the A jacobian is shown here: 
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with the following terms: 

421 = o/3 — u 2 
A22 = (2 - (3)u 

441 = 20au — yEu + /?& — vu v” — tin” 2 

3?1 2 4- f> 2 — 

442 = 7 E- /3Jfc + u ” 2 
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These matrices show clearly the coupling between the Navier-Stokes equations and the 
transport equations for the normal components of the Reynolds stress tensor, just because 
of the introduction of the turbulent kinetic energy in the global energy budget, whereas 
there is no apparent coupling with the shear stress and the dissipation rate equation. In 
fact, these equations are related to the previous through the source terms only. 

Ill- 2-3. Treatment of non-conservative equations 

To treat implicitely the source terms, various techniques are available. Recall first the 
general implicit approximation: 

(I + A + A - A t.C) 6U n+1 = A U n ( III - 25) 

ox ay 


25 


The simplest way, which is somehow trivial, is to apply a first approximate factoriza- 
tion, without considering the formal content of the source terms. Then it comes: 

{I+At ^ir +At i£- At - C)={I+At ^ 7 + - a<.c) + o< a* 2 ) (/// - 2 s> 

The C matrix can be considered as diagonal whatever the formal content is, i.e. 


C = 


ro o o o 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 

0 0 0 0 
0 0 0 0 

0 0 0 0 
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Before doing the work on the space operator, the inversion of the diagonal source term 
matrix is straightforward: 


(7 + At^ + Af^) 6U n+1 = AU n . (7 + A*|C|) 


-l 


(777 - 28) 


Therefore, the explicit increment is modified first by the source term contribution, 
before being updated by the space derivative operator(s), either with an approximate 
factorization or a relaxation technique. 

A slightly different method avoids the factorization for the source contribution. Then 
the source terms are grouped with the transverse advection operator [20], [21]. In that case, 
the same eigenvalue is used, which is the the maximum value among all equations to be 
solved. 

Unfort unately, the use of these blind forms, without accounting for the formal content 
of the source terms does not guarantee the stability. Therefore, it has been found necessary 
to develop more exact forms of the jacobian matrix. Although various developments are 
possible, we will develop here a typical form which has been prooved very efficient, as far 
as stability is concerned. 

Consider the set made only with the turbulence transport equations. The convective 
and diffusive parts axe supposed already solved with the Reynolds averaged Navier-Stokes 
equations. Then we have only to work on: 


dU 

dt 


= H 


( III - 29) 
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where the two vectors U and H are now: 


U = 


p k 
pt 


H = 


H k 

H f 


{III - 30) 


An implicit approximation of equation (III-29) is: 

U n+1 = U n + At. H n+1 


{III -31) 


in which H n+1 is evaluated at time {n + 1). This can be achieved by a first order serie 
expansion: 

with SU = U n+1 - U n 


ff n+l =H n + <W 6U 

oU 


Then the implicit approximation can be rewritten as: 


{I - At SU = At H n {III - 32) 

v ou 

The task is to evaluate properly the jacobian matrix. Let first discriminate between 
positive and negative source terms. An rather elementary stability analysis on equation 
(193) shows that stability cannot be obtained with an implicit approximation when the 
source term is positive. The same is true for an explicit approximation with negative 
source terms. Therefore we keep only in the implicit approximation the ’’good” terms, 
which are negative. All the permanently positive contributions are treated exclusively in 
the explicit part of the scheme. 
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IV. CONCLUSIONS 


In the framework of this report, the following tasks have been accomplished: 

* Implementation of variety of transport equation turbulence models in a versatile 
boundary layer code. These models range from the algebraic Baldwin- Lomax to the 
full second order closure, derived from the LRR approach. 

* Implementation of a full second order closure in an implicit solver (MacCormack 
scheme complemented with flux vector splitting and line Gauss-Seidel relaxation 
method. 

These various turbulence model implementations have been applied to a wide range 
of compressible flows in two dimensions. 

The second order closure have been shown to account implicitely for complex turbu- 
lence effects, like strong anisotropy variations or curvature effects. Nevertheless, experience 
of the authors have shown that its use for routine computation is still limited by stiffness 
of numerics and computational costs. 
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Part II. 


C. Dingus and W. Kollmann, MAME Dept., UCD, Davis, CA. 95616 


Objective. 

The objective of the present part of the project was the development of a complete second 
order closure for wall bounded flows including all components of the dissipation rate tensor 
and a numerical solution procedure for the resulting system of equations. The main topics 
of the present grant were the closure of the pressure correlations and the viscous destruction 
terms in the dissipation rate equations and the numerical solution scheme based on a block- 
tridiagonal solver for the nine equations required for the prediction of plane or axisymmetric 
flows. 
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1.0 Full second order closure for wall bounded flows. 


1.1 Introduction. 

A full second order closure for wall bounded shear flows is developed, which includes the 
Reynolds stress equations and the equations for all relevant components of the dissipation 
rate tensor. Incompressible and compressible plane flows are considered, but this report is 
only concerned with incompressible flows. 

There are several reasons for the development of a complete second order closure. It can 
be shown that the anisotropy of stress and dissipation rate tensors approaches the same limit 
at the wall, but the derivative of the anisotropy of the dissipation rate is twice the derivative of 
the stress anisotropy at the wall. Another reason is the possibiity of constructing apppropriate 
time scales in the near wall region. The standard second order closure incorporates the 
transport equation for the trace of the dissipation rate tensor and relates the components of 
the tensor via local relations to the trace. The time scale for the destruction of the trace is 
usually modelled using the time scale 



e 

k 


where the modified dissipation rate i is defined by 



with y denoting the coordinate normal to the wall. This is an acceptable model in the region 
close to the wall if and only if the dissipation rate is a nondecreasing function of the distance 
from the wall, because the kinetic energy is of order 0{k ) = t/ 2 near the wall which implies 
that the second term in the modified dissipation rate is constant. Direct simulations of 
boundary layers and channel flows, however, have shown (see Mansour et ah, 1988) that the 
dissipation rate is a rapidly decreasing function of the wall distance in the viscous sublayer 
and the destruction model using the time scale tq becomes thus a production term in the near 
wall region. This is clearly a violation of realizability for the destruction model. It follows 
that this type of closure model does not represent properly the production of dissipation rate 
near the wall and the sign reversal of the viscous destruction model must make up for this 
deficit in production. The motivation for the modified time scale is the fact that the time 
scale 

_ * 

e 

goes to zero as the wall is approached. It will be shown below that the full dissipation rate 
tensor allows a realizable and tensorially invariant construction of a time scale or a time scale 
tensor that reaches a finite and nonzero limit value at the wall. 
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A second point that sets the near wall region apart from the high Reynolds number 
regime of the boundary layer is the growth of the pressure corelations with distance from the 
wall. It will be shown below that the usual split of the pressure correlations into pressure 
transport and pressure rate of strain correlations is not appropriate near the wall, because the 
split correlations grow with different rates and the Taylor series for the original correlation 
involving the fluctuating pressure gradient can be expressed locally up to second order in 
terms of velocity correlations. 


1.2 Exact equations for the dissipation rate tensor. 

Incompressible flows are considered first. Standard manipulations lead to the transport 
equations for the dissipation rate tensor in a Cartesian coordinate system defined by 

e a p = 2ud- t v' a d-,v , 0 (1) 

with expectation denoted by {e Q p). The equations can be given in the form 

( d t + (i; 7 )d 7 )(e Q/ ?) = d 7 [t 'd 7 (e Q /j) - ( u 7 e a/?)] + Sip + S 2 a p + Sip + S*p + Sip - D a p (2) 


The various source terms are defined as follows. 

Sip = —(ep^dffoa) — { cq 7 )9 7 ( u ^} (3) 

Slp = -(el 6 p)(d y (vs)+ds(v y )) (4) 

Sip = -(e^d y v' Q ) - (e a ^v'p) (5) 

SU = -j((dl y p'd^) + (d^p'dX)) («> 

Sip = -2u((v' 6 d y v'p)dl s {v a ) + (vsd^v' a )d 2 yS {vp)) (7) 

D a p = Au 2 {d 2 h v' Q d 2 6l v , p) (8) 

Furthermore is the fourth order dissipation tensor defined as 

*tp = 2vd 6 i >' a d y v'p ( 9 ) 


The properties of the source/sink terms on the right hand side need to be established before 
closure expressions can be constructed. 

1.2.1 Equivalent forms of the pressure correlations. 
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The pressure correlations can be given in several equivalent forms. It is instructive to 
split them analogous to the pressure-rate of strain correlations in the stress equations. It 
follows from (6) that 


S^p = -—{d Q (diP'd- t v' l3 ) + dp (d^p'd y v' a )} + B a p 

where the non-gradient part is defined by 

B a p = ^-(d^'d^dcv'p + d fi v' a )) 

The non-gradient part B a p resembles the pressure-rate of strain correlation and shares with 
it the property 

B qq = 0 

It follows that B a p redistributes intensity among the components of the dissipation tensor 
and leaves its trace unaffected. Further splitting of B a p leads to 

B a p = —d 1 {d 1 p'{d Q v , 0 + dpv ' a )) (d* y p' (d a v'p + dpv' a )) 

P P 

The non-diffusive part of B a p contains the Laplacian of the pressure fluctuations which is 
governed by 

-~d*p' = dfiv' a d a V 0 + 2 d a v' 0 dp(v a ) - dlp(v' Q v'p) 

P 

This equation has an important consequence: The non-diffusive part of B a p can be repre- 
sented locally in terms of velocity fluctuations. We get 

B a p = —dfid-jp'idav'p + dpv' a ))+ 

P 

2u (dyVfidsv'^daV'p + dpv'a)) + + dpV Q ))ds(v 7 ) 

We conclude that the non-diffusive and non-gradient part of the pressure correlations does 
not contain a direct influence of the wall and the wall effects can be represented as gradients 
and divergence of a flux. This property very important for the modelling effort. Inspection 
of the local part of B a p shows that it has the structure of the primary production term (5). 
Recasting this part in terms of vorticity and strain rate defined by 


and 


s aP = ~^{da v P + Bp v a) 
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leads to 

2u(d~ 1 v'fid6v'^(d a v'i3 + dpv'a)) = 4l / { s Q p S y6 S ts) ~ 
whereas the primary production terms appear as 

Sip = —^({ s 'a-r s 'iS s 'sp) + € iQr t {s' l 3is' 6l uj' 11 ) + e 1 p^(s' QS s' s ^u>' TI ) + e 1QU ,esp v {s yS u ^ J)) 

The last term can be recast in terms of the Kronecker delta using the tensor relation 


/ b-,6 

e^ocw^-S pi) — (let I Safi 

\tiu>6 


fiyP 

&ap 
6u >p 


^OIJ 


leading to 


sip = -^({s'ays'yss'tp) + e-ra V 

+ ( s '-,p Lj 'a u> 7) + {s'-ia^'p^'-i) - ( 


($ Pi^S-f^ij) “b €-)PTl{ S ot6 S Sy^i]) 
S' a pUj'yJy) ~ 8ap{s'y^ V'yU'j) 


It is apparent that no complete cancellation of triple correlations takes place. The pressure 
correlations appear now as a combination of the divergence of a flux and sources. 


Sit = ft**,. + Qlt 


where the flux is defined by 




{eujSrjFutyaidyP drjVp) -f- (^~tP ^V V a)') 

P 


Kps - -—({dsp'idav'p + dpv' a )) - Sasidyp'dyv'p) - f>ps{dyp'dyv' a )) 
The trace of the flux F^ s is not zero but given by 

C1-7WW] 


The redistributive sources are given by 

Q% = 2,/(3 1 i ’'Av'psyt + deO) + ^(SA(9av't + dtv' a ))d,(vj 

where the trace of Q p a p is zero. It is noteworthy that one of the components of (d-^p d-qV a ) 
can be expressed in terms of a component of the dissipation tensor 

4- (chp 1 chv'J - vd-ifai) + 0(y) 

P 
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according to the power series for the pressure given below (19). For 7^2 the gradient of the 
wall pressure can be expressed in terms of velocity using the momentum balances. In fact 

- p d -*p' = v8 iA 

holds at the wall y = 0.0. It follows that the wall pressure does not exert a direct influence on 
the dissipation rates if the expansion is carried out to second order. It is clear that modelling 
can be based on the properties of the flux and the redistributive source terms. However, the 
growth rates for the terms in the different formulations of the pressure correlations decide in 
the end their usefulness. This will be investigated in the following chapter. 


1.2.2 Taylor series expansions for the near wall region. 

The near wall region can be analyzed with Taylor series. The coordinate system is 
assumed to be located at the wall and x 2 is the direction of the wall normal pointing into 


the flow field. It is convenient to rename the coordinates and variables as follows: xj = x, 
X 2 = y, X 3 = z and 17 = u, v 2 = v, t > 3 = w. The velocity components can be expanded with 
respect to the wall normal y 

u(x, y, z,t) = a 0 + my + a 2 y 2 + a 3 y 3 + 0(y 4 ) (10) 

v(x,y,z,t) = 60 + hy + hy 2 + b 3 y 3 + 0(y 4 ) (11) 

w(x,y,z,t) = c 0 + cnj + c 2 y 2 + c 3 y 3 + 0(y 4 ) (12) 

where the coefficients are stochastic functions of x,z,t but not y. They are defined by 


The noslip condition at the wall 

1 d 3 u 
1 d J w 

c > (x, 2 ,() = i , a ^(°) 

implies that 



a 

0 

II 

0 

II 

II 

0 

(13) 

holds and mass balance 


0 

II 

JS 

a 

(14) 

leads to 


d y v 0 = bi = 0 

(15) 
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and 



d x a n + d z c n = -(n + l)6„+i 


for n = 1,2, • • 

•. The Taylor series for a fixed wall without suction or blowing are 

therefore 

given by 

u(x,y,z,t) = a 3 y + « 2 y 2 + a 3 y 3 + 0(y 4 ) 

(16) 


v(x,y,z,t ) = b 2 y 2 + b 3 y 3 + 0(y 4 ) 

(17) 


w(x,y, z,t ) = cry + c 2 y 2 + c 3 y 3 + 0(y 4 ) 

(18) 

The expansion for the pressure can be given as follows 



p(x,y,z,t ) = po + 2 pb 2 y + 3 pb 3 y 2 + p 3 y 3 + 0(y 4 ) 

(19) 


where /z denotes the dynamic viscosity and 

Pj{x,z,t) = 

and the momentum balance normal to the wall was applied to a point at the wall. The 
expansion for the pressure gradient can be shown to be 

(a 2 \ /di& 2 \ (fidxbA 

d-yp = 2p I 62 ] +2 py j 3&3 ] + 3y 2 I p 3 1 +0(y 3 ) 

\ c 2 / \d 3 b 2 ) \ pd 3 b 3 J 

which shows that the terms up to first order are proportional to viscosity. The components 


of the Reynolds stress tensor appear in expanded form as 

(it 2 ) = y 2 (a 2 ) 4- 2y 3 (aia 2 ) + y 4 (2(ai<x 3 ) + (a\)) + 0(y 5 ) (20) 

(v 2 ) = y i {bl) + 2y 5 (b 2 h) + 0(y 6 ) (21) 

( 10 2 ) = y 2 (c 2 ) + 2 y 3 (c lC2 ) + y 4 (2( Cl c 3 ) + (c 2 )) + 0(y 5 ) (22) 

(uv) = y 3 (ai6 2 ) + y 4 ((ai& 3 ) + (a 2 & 2 )) + 0(y 5 ) (23) 

The components of the dissipation rate tensor vary near the wall according to 

(e n ) = 2is{(a 2 } + 4y(«ia 2 ) 

+y 2 (6(a ia3 ) + 4(a 2 ) + ((d^) 2 ) + ((3,a,) 2 )) + 0(y 3 )} (24) 

(c 22 ) = 2v{\y 2 {b\) + 12y 3 (6 2 6 3 ) + 0(y 4 )} (25) 

(£33) = 2 ^{(c 2 ) + 4 y(cic 2 ) 
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+y 2 (6(c!c 3 ) + 4{c 2 2 ) + ((«W> + (&ci)) + 0(y 3 )} (26) 

(^ 12 ) = 2^{2i/{ai6 2 ) + y 2 (4(a 2 6 2 ) + 3{ai& 3 )) + 0(y )} (27) 

It is apparent from these series that different components grow differently near the wall. It 
follows that the boundary values for the dissipation rate tensor are given by 

<e n )(0) = 2i/(a?) (28) 

(e 22 )(0) = 0.0 (29) 

< c 3 3 >( 0 ) = 2 u{c\) (30) 

(e 12 )(0) = 0.0 (31) 

It can be shown that the anisotropy of the Reynolds stress tensor is equal to the anisotropy 


of the dissipation rate tensor at the wall and that the normal derivative at the wall of the 
anisotropy tensor of the dissipation rate is twice the normal derivative for the stress ten- 
sor. We consider now the near wall variation of the individual terms in the dissipation rate 
equations. 


1.2. 2.1 Viscous Diffusion. 

The dominant term in viscous diffusion is the normal derivative given by 

dy(vd y {en)) = 4i/ 2 (6(aia 3 ) + 4 (a]) + {(d^) 2 ) + (( d *aj) 2 )) + 0(y) (32) 

d y (ud y (e 22 )) = 16i/ 2 (6 2 ) + 0(y ) (33) 

dy( u dy{ e w)) ~ 4^ 2 (6(cic 3 ) -I- 4{c 2 ) + {(dx c i) 2 ) + {(d* c i) 2 )) + 0(y) (34) 

9y{t / dy(e 12 )) = 4i/ 2 (4(a 2 6 2 ) + 3(ai6 3 )) + 0(y) (35) 

and they emerge as terms of order unity near the wall for all components. 


1.2. 2. 2 Turbulent diffusion. 

Turbulent diffusion of the dissipation rate component e a p is defined by 

d y Ff = -d y {v' y e a0 ) (36) 
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For the case of boundary layer type flows only the flux normal to the wall is relevant. The 
series expansions lead then to the following expressions for the components of the dissipation 
rate tensor 


d y F" =2y(b 2 e u (0)) + O(y 2 ) 

(37) 

d y Fy 2 = 4j/ 3 (i>2f22(0)) + 0(y 4 ) 

(38) 

d,F™ = 2y(bi' 33 (0)) + 0(v 2 ) 

(39) 

S,e; 2 = 3 v 2 ( 62 a,ei 2 ( 0 ))+O( ! / 3 ) 

(40) 


The expansions show that turbulent diffusion is not of leading order near the wall. 


1.2. 2.3 Secondary production S^. 

The interaction of the mean rate of strain and the dissipation rate tensor acts as production 
for the dissipation rates in the same fashion as the Reynolds stress and mean strain rates 
for the stresses. There is however a fundamental difference between this type of production 
for dissipation rates and stresses: It is of leading order for the stress balance but of second 
order for the dissipation rates for high Reynolds number flows. The situation near the wall 
is entirely different. The components of turn out to grow with wall distance as follows 


Sh = -2yd y (u)(0)d y (enm + O(y 2 ) 

(41) 

SJa = -16^ f {t;)(0)(^)+O(y 4 ) 

(42) 

w 

II 

O 

O 

(43) 

5 1 1 2 = -8 I /y 2 5 J/ (u)(0)(6^)+O(y 3 ) 

(44) 


The secondary production terms are not of leading order near the wall, but closed. Hence, 
they need not be neglected. 


1. 2.2.4 Secondary production S 2 a p. 

The series expansions lead to the following results 

Sji = -y{a,(u)(0)a,(e 11 )(0) + 2Sj # (^)(0)(€„)(0)) +0(!/ 2 ) (45) 

Sj 2 = -2vy 3 {2d,(u)(0)d,(b 2 ) + 8aJ„{.>)(0)<& 2 )} + 0(y A ) (46) 

Sf, = -vW<*>(0)3«(«33)(0) + 2a 3 ,(t>)(0)(€„)(0)} + 0(v 2 ) (47) 
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S? 2 = -2i'!/ ! {2a s {u)(0)((6 2 a i o,) + (a,d,h)) + 4d^(r}(0)(a 1 b,)} + 0(y 3 ) (48) 

It is clear from these expansions that the source terms S are of second or higher order for 
boundary layer flows. 

1.2. 2. 5 Primary production. 

The primary production or vortex stretching terms are the dominant production terms in 


high Re-number flows. For the near wall region they appear in expanded form as 

Sn = -4 v[y{\dM) + + 0(y 2 )) (49) 

S| 2 = — 4i'{y 3 (10{6 2 ) + a,(a,6|) + 0(y 3 )) (50) 

= -4*{y(ia,(a,c?) + 3(c]h)) + 0(y 2 )} (51) 

Si 2 = ~'2v{y 2 {d x {b 2 a\) + 8(ai b\) + 2 (c x b 2 d z a x ) + {a x c x d z b 2 )) + 0(y 3 )} (52) 


The primary production is not of leading order near the wall, but grows with the same order 
as the secondary production term 5^ with wall distance. 

1.2. 2.6 Viscous destruction. 

Viscous effects can destroy the rate of dissipation and this process is contained in D a p. The 
series analysis shows that the components of D a 0 are near the wall given by 


D n =4y 2 {4(a|) + 2(O,a 1 ) 2 ) + 2((a,a 1 ) 2 )}+0(y) (53) 

D 22 = 4i' 2 {4(6 2 ) + 24y(6 2 6 3 )} + 0( y 2 ) (54) 

£> m = 4i/ 2 {4(c 2 } + 2((a,c, ) 2 ) + 2((3,c, ) 2 )) + 0(y) (55) 

Du = 4i/ 2 (4(a 2 4 2 ) + y(4(d r ai3 r 5 2 ) + 4{d : a 1 d : & 2 )+ 

12(a 2 6 3 ) + 12(<ij6 2 )))+0(y 2 ) (56) 

All components of D a p turn out to be of leading order. 

1.2. 2. 7 Pressure correlations. 
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The Taylor series for the pressure (19) contains viscous terms which are due to the momentum 
balances. Differentiation of (19) leads to expansions which contain viscous contributions in 
lowest order. The components of the pressure correlations S* 0 appear as follows 


Sf x = 4i/ 2 {6(ai6 3 ) + y{2{d x a l d x a 2 ) + 2(d s a 1 d z a 2 ) + 4(a 2 d x b 2 ) + 6{ai<9 r 6 3 )) + 0(y 2 )} (57) 

S\ 2 = 4i/ 2 {12 y{b 2 b 3 ) + 0(y 2 )} (58) 


5, 3 = 4 i / 2 {6( ci 6 3 ) + y{2(d x c\d x a 2 ) + 2(d z c l d z b 2 ) + 4{c 2 d x c 2 ) + 6(c { d x h)) + 0(y 2 )} (59) 


6 

s ; 2 = 2i/ 2 {6{a,ij)+!/(2(a t aia t 6 2 ) + 2(a ; a,a,ii2) + 12(a 2 i> 3 ) + 4 {M«M + — (a t p 3 )) + 0(y'‘)} 

(6°) 

The series expansions of the pressure correlations lead to several important conclusion: The 
effect of the pressure correlations in lowest order is local in terms of velocity correlations. No 
Poisson integral appears in lowest order since the wall pressure does not appear in the lowest 
order terms. The split of the pressure correlations obtained in chapter 2.1 leads to a flux 
such that the corresponding source term is strictly redistributive and local in terms of 
velocity fluctuations. The pressure flux is local in velocity in expanded form up to second 
order. Note that there is a clear advantage for not splitting the pressure correlations in the 
Reynolds stress equations near the wall since the split terms (rate of strain correlation and 
pressure transport) grow with different rates in the viscous sublayer. This is not the case for 
the pressure correlations in the dissipation rate equations. The components of the pressure 
flux emerge for the case of a flat plate boundary layer as 

4- u 

f'u = j&p'bv ;> 

F? JS = — {(hp'dWt) + (8iP'M)} 

P 

At/ 

p 

F r m = j {Wdnu - (d 3 p'd*v\)} 

which can be analyzed with the aid of Taylor series. We get the following estimates for the 
divergence of this flux near the wall 

diFf n =Sv 2 (b 2 d x a x ) + 0(y) 
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d 2 F > 22 = 0(y) 

d 2 Fl v = 8 v 2 {b 2 d lCl } + 0 {y) 
chFf 22 = -4^ 2 {c 2 c) : c 1 ) + 0(y) 

It follows that 

diFf u + d 2 F ^ 2 = -d 2 (ud 2 (e 22 )) + O(y) 

holds. This proves that the pressure flux terms are of leading order near the wall because 
viscous diffusion is of leading order in this region. 


1.2. 2. 8 Near wall production. 

The production terms proportional to the curvature of the mean velocity profile are respon- 
sible for additional production in the near wall region. The series expansions lead to 

Sj, =^{y{a;>a^(u)(0)+y 2 ((o;)^„(“}(0) + 3(<i 1 a 2 )a^( u )(0) + 2(a 2 i 2 )9j,(u)(0)) + 0(^)> 

Sf 2 = 0(v 5 ) (62) 

Sf 3 = 0 (63) 

S ? 2 = K» J {(«?>fl{„(«)(0) + (ai6j)^„<«)(0))) +0(y 3 ) (64) 

The near wall production emerges as second order effect near the wall if the boundary layer 
assumptions are satisfied. The component (e n ) receives all the energy in lowest order. 


1. 2.2.9 Transport equations in lowest order. 

The series expansions for the source and diffusive terms allow the set up of the transport 
equations for the components of the dissipation rate tensor in lowest order. It turns out that 
all equations are of the same zeroth order. 


<9({en) = d 7 (Vd 7 {eii)) - D n + Sn 

(65) 

0 = d- 1 (v / d- 1 (e 22 )) — D 22 

( 66 ) 

dt(e 33 ) = 5 7 (^5 7 (f33)) — D 33 + S 33 

(67) 

0 = d- y (i , d 1 (en)) — D \ 2 + Sf 2 

( 68 ) 


These equations are valid near the wall provided none of the surviving correlations vanishes. 
For steady flows the following equations hold then at the wall 

D n -5?, =d y (ud y (e n )) 

D 22 = d 1 (ud~ l (e 22 )) 

I>33 - S 3 4 3 = 07(^33)) 

D u - SU = d^ud^en)) 

which can be combined with the limit relations for the stress balances to establish constraints 
for the modelled terms at the wall. 
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1.3 Closure model for the dissipation rate equations. 


The series expansions for the near wall properties of the dissipation rate equations can be 
used to analyze and to modify closure expressions. In several cases no model expressions 
exist and new models will be developed and analyzed. 


1.3.1 Time scales. 


Several time scales can be constructed for the near wall region with the aid of the dissipation 
rate tensor. First we note that a scale dependent on the wall normal vector can be obtained 
in the form 

r = (69) 

l~l ff ^ 6 $ f 

If n y = S-,2 it follows that this time scale is given by 

k 


T = 


( e 22) 


and the series expansions show that both numerator and denominator depend quadratically 
on wall distance. The wall limit is in fact a nonzero value given by 


, . _ ( [dyUodyUQ ) + {dyWpdyWo) 
T ~ ^{0 2 yy V 0 dl y VQ) 


(70) 


Hence, a time scale with a nonzero limit at the fixed wall was constructed. This time scale 
avoids the problem associated with the modified (also called homogeneos) dissipation rate 

e = e - 2u{d y \fkf 

which may change sign in the flow field. The inverse of another time scale with tensorial 
character using the dissipation rate tensor can be set up as follows 


T al = l {e 6 p) 


( 71 ) 


where (v l y v' 6 )~ 1 denotes the inverse of the Reynolds stress tensor. Conversely is a time scale 


given by 

Tap = ( V 6 V 'p) 

where {e a /}) -1 denotes now the inverse of the dissipation rate tensor. 


( 72 ) 


1.3.2 Turbulent diffusion. 
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The gradient flux model for can be given in the form 



k 

e 


(v\v' 5 )d 8 {e a o) 


where the kinetic energy is denoted by 


( 73 ) 


k(x,t) = ^{v’ a v' Q ) 

and the dissipation rate e by 

e(x,t) = ^{e oa ) 

The constant c. s has values in the range 0.15 — 0.18. The near wall properties of this model 
follow from the series expansions (10) to (27) as 

d y F"=0(y 5 ) 

whereas the exact term has a first order variation with respect to the wall distance according 
to (37). Similar discrepancies are observed for the other components. It is clear that model 
expressions developed for the high Re-number regime will not represent the near wall region 
properly. The present model (73) implies that near the wall turbulent diffusion is essentially 
neglected in comparison to the exact term. Inspection of the model (73) shows that there are 
two reasons for its failure near the wall: The time scale * approaches zero at the wall and 
the diffusivity is solely determined by the normal stress component ((V) 2 ) which varies as y 4 
near the wall. The situation can be improved if a composite time scale that approaches the 
scale defined by (69) near the wall is used. The modified closure model is then given by 


F a ^ = 

7 



k 




(74) 


The factor 2/3 results from the requirement that the high Re- number limit of the time scale 
must agree with k/e. 

A different model that satisfies all growth estimates and has the correct tensorial and 
dimensional properties can be constructed if the dependence of the turbulent flux on the wall 
parameters n and i? e < is taken into account and the near wall model is combined with the 
high i? e -number model. The near wall model is given by 

-{v' y e a ^)=c'J w ■2(Ret)[n f) n < (v' v v' ( .)}^n y (e Ql 3) 

It is straightforward to check that it has the same growth rate and the same tensorial prop- 
erties as the exact term. The model represents turbulent transport towards the wall and 
has the form of a convective term. There exist several high Reynolds number models for the 
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turbulent flux of dissipation rate (Hanjalic and Launder, 1972, Lumley, 1978). The present 
model is an analogue of the flux model for the stresses. It is given by 



d 

dxs 


(cq/j) 


The combined closure model is the set up as follows 

Jc d 

-(v' i e Q 0 )=c'J w 2 (Ret)[nr l n (: {v l rj V l < .)}h-l y (e Q 0 ) + c e (l - fw 2 (Ret))-(v\v' 6 ) — ( 6 ali ) 


where f W 2 (Ret) denotes a function of the local Reynolds number such that f W 2 goes to unity as 
the wall is approached and to zero in the turbulent zone. Furthermore, The function fw 2 (Ret) 
should be nonnegative and it should not modify the dependence on the wall distance for the 
near wall model. It follows that f w 2 must be an exponential function of the Reynolds number 
given by 

fw2{Ret) = 

Ke2 

where R e 2 is a constant measuring the range of influence for the near wall model. This 
function of the Reynolds number has the well known property that all its derivatives at zero 
vanish. Hence, it does not modify the growth rate of the near wall model. 


1.3.3 Secondary production S 2 a0 . 

The properties of Sl 0 in the high and low Re-number limits will be considered first. Local 
isotropy requires that 

„ lim se 

Re— >00 y 

holds. It follows that the secondary production terms S^ 0 are for high Re-numbers given by 

lim S 2 a0 = -^ a /t A A) (75) 

Re— ^00 y 

The divergence of the mean velocity is zero for incompressible flows and it follows that the 
high Re-number form of the secondary production terms can be neglected. The near wall 
variation of the secondary production term 5^ according to (45) can be regarded as 

Sf, = -(u)3,<«i,) - 2(«„)3,(»> + 0(y 2 ) 

and it follows that this term is at best of the same order of magnitude as the mean convection 
term. Similar relations hold for the other components and we conclude that the secondary 
production terms S^ 0 can be neglected. 

1.3.4 Primary production and viscous destruction. 
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The closure model for these terms is fundamentally different in the high and the low Reynolds 
number limits. Both limits need to be considered and the corresponding closure models must 
be merged to cover the range of Reynolds numbers from zero to infinity. 

High Reynolds number limit. 

The order of magnitude estimates for primary production and viscous destruction at high 
Re-numbers shows that they are of leading order, but their difference is of the same order 
as the secondary production term S l Qj3 . It follows that they should not be treated separately 
but their difference should be modelled as function of the available information. The model 
consistent with second order moments is in general given by 

Sip - D a p=-<lf Q p({€ Sul ), {vev'JydsM) ( 76 ) 

The dimensionless and symmetric tensor should represent both production due to the 
interaction of vorticity and strain rates and the destruction due to viscous effects. If we 
impose the condition of local isotropy on this model we get the following variant 

dy(v$)) — { v s v <jj )> ^ i v u) ) 

where the first part represents the productive and the second the destructive contribution to 
the model for the difference of 5^ and D Q @. The model can be set up to be consistent with 
the standard expressions for the trace equation (see Launder, Reece and Rodi, 1976) 


* aj3 — ^ 


d-t{ v s) - c (2 


i^p) 

e 


and the time scale is given by 


The closure model 


T = 


k 

e 


~ D Q f)— — C t i ^ap ^ { V f V s)^‘l( v ^) c (2{^ap) ^ ( 77 ) 

emerges. It is, however, not applicable to wall bounded flows since the destructive part of 
this model becomes singular as the wall is approached. This deficiency can be corrected 
either by defining a time scale that does not go to zero at the wall or by merging low and 
high Reynolds number closures with a Reynolds number dependent function such that the 
singularity is removed. Finally, we note that this model is not necessarily positive definite 
because the productive part d /(1 * is not positive definite, a property shared with the exact 
term. Refined closure models can be constructed using tensorial time scales introduced in 
chapter 1.3.1. 

Low Reynolds number limit. 
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The dependence on the wall distance is the deciding property as the wall is approached. We 
recall that the primary production terms decay with wall distance, according to chapter 2.2.5, 
as follows 

0(S 3 U ) = 0(5 3 3 3 ) = y 
0(S 2 ’ 2 ) = y 3 
0(S? 2 ) = y 2 

whereas the viscous destruction terms are all of leading order 


0{D Qt )) = 1 


according to chapter 2.2.6. It follows that they must be modelled separately in the near wall 
region. The first step in the construction of the low Reynolds number version is the analysis 
of the near wall properties of the high Re-number model. It follows from (77) that 


0(^ {l) ) = y 


which is at variance with the detailed decay laws for 5^ given above for a = f3 = 2 and 
a = 1,/? = 2. However, since the primary production is not of leading order near the wall 
it would be acceptable. A more serious problem arises in the destructive part ^ a p. If we 
construct a time scale such as (69) for the near wall region we avoid the singularity, but it is 
not possible to satisfy the decay laws for the destructive terms D a p. It follows that the time 
scale in any closure of the form given by cannot be a scalar but has to be a tensor of 
rank two (or higher) with positive eigenvalues. It is not difficult to construct a closure model 
such that the decay law is satisfied. For instance, the model 


Daf) — 


c e2 r 


6 


deaf dtfff 
dXfidXj! dXfjdXq 


with (69) as time scale possesses the correct tensor properties and the correct decay law as 
the wall is approached. However, it is unacceptable as closure model because it produces 
in regions where the second derivatives are all positive a second order pde with negative 
diffusivity. The initial/boundary value problem for such an equation is not well posed and 
the numerical solution futile. Furthermore, this model would not be realizable, because the 
limit (e a p) 0 does not imply that the second derivatives of {tap) go to zero. It follows that 
the closure model for the viscous destruction of the dissipation rate components must be of 
the form 

D a(i =F Ql i{{t v ),{v\v' 1] ),d (> {tv ] ),d ( ' (vyj^Ret) 
where R et denotes the turbulent Reynolds number defined by 

k 2 

Ret = - ( 78 ) 

eis 
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The example cited above for the destruction of dissipation rate showed that derivatives lower 
that second should be used. Several closure models can be set up that are tensorially consis- 
tent, have the correct decay law near the wall and avoid the stability problems incurred by 
negative diffusivities. We note two of them, first 


D a f}= — C ( 2 


r 


C e 3 T TlfTlfi 


dx y dxs 


(79) 


and 


D a jj — 



c e3 u€7i y ns 



( { °i) \ d r 
6 } dxs y 


e 


(80) 


as second variant. Note that n Q denotes the unit normal vector of the wall pointing into 
the flow field. The time scale r is given by (69) to avoid a singularity at the wall. The two 
models look very similar, but inspection of the time scale (69) shows at once that the first 
model is most likely unstable. Suppose that kinetic energy k and dissipation rate component 
n 1 ns{e 1 s) = ( 622 ) are in equilibrium, then assume that this equilibrium is disturbed, say the 
kinetic energy is reduced by a small amount. It follows from the fact that the second part in 
the first model is quadratically proportional to the kinetic energy that the rate of destruction 
of the dissipation rate components is decreased by the disturbance and consequently are the 
normal stresses further decreased and the equilibrium state is not recovered. Hence is this 
model unstable. It follows that the second model is the preferred one. 


Merger of low and high Re-number models. 

The closure models for the destruction term (79) or (80) valid for the limit R e -» 0 and 
the destructive part of (77) valid for the high Re number limit must be merged together to 
produce a model valid throughout the flow field. Suppose both limit expressions have decay 
laws near the wall that do not need change via the function weighting them according to the 
local Reynolds number (or any other function propertional to the distance to the wall). The 
we need a weight function that does not change the dependence on wall distance near the 
wall. This implies that we must find a function which has zero derivatives at zero. It is well 
known that the exponential function 


/(y) = ex p(-i) 

y 

has the desired property, in fact, it is infinitely often differentiable but not analytic at zero 
(its Taylor series is identically zero at y = 0). Hence, we can establish a low Reynolds number 
function po 

f(R e ) = exp[-(-^-) 2 ] 

which is zero at zero Re-number and unity at infinite Re-number. The constant R t determines 
the range of Reynolds numbers for which the function is close to zero. Other functions 
have been proposed that are proportional to some power of the wall distance and change 
therefore the decay law. Several models have been proposed for the functions f \ , multiplying 
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the production of dissipation rate, and / 2 , multiplying the destruction of dissipation rate. 
Vandromme et al. (1983) suggested 

fi = 1.0 


and _ 

Tie 

f 2 = 1.0- 0.22 expM^) 2 } 

based on the original model of Hanjalic and Launder (1976) and obtained good agreement 
with measurements in flat plate boundary layers. Recent developments surveyed by Launder 
(1989) use the second and third invariants of the Reynolds stress tensor to represent the wall 
influence on production and destruction processes. 


1.3.5 Pressure correlations. 


It was shown in chapter 1.2 that the pressure correlations can be split into transport and 
source terms such that the source terms are strictly redistributive and have no effect on 
the trace of the dissipation rate tensor in analogy to the pressure-strain correlations for the 
Reynolds stress tensor. It follows that they must redistribute intensity among the components 
of the dissipation tensor. Kolmogorov’s hypothesis of local isotropy requires that the dissi- 
pation rate tensor approaches its isotropic form as the Reynolds number approaches infinity. 
It follows that a return to isotropy model given by 

7« e «o> - !«■*') < 79 ) 

would satisfy this condition. The open question is the time scale. The scale 

- t 

6 

is the obvious choice for the high Re- number limit, but the model becomes incorrect as the 
wall is approached because r -1 goes to infinity with j/ -2 . Modification of the time scale 
according to (69) solves this problem and 


( 80 ) 

emerges as nonsingular variant. The transport part of the pressure correlation can be mod- 
elled in terms of the viscous diffusion terms because the wall limits indicate this form. The 
model 

^'ai3-f = ^Cctrjui^ ))^ui ~f" £ Or } (81) 

where the tangential unit vector is defined by 


= M) 

“ \M(6)\ 


(82) 
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and the binormal by 


fra = €a fifty 


satisfies the near wall properties for the components F[ 12 and Ff 32 and neglects the compo- 
nents F% 22 and F[ 22 . The effect of the diffusive part of the pressure correlations is therefore 
inhibition of the viscous diffusion near the wall for the diagonal components corresponding 
to motion parallel to the wall. The effect on shear component and the diagonal component 
corresponding to the motion normal to the wall are neglected. 


1.3.6 Near wall production. 

The near wall production requires in general flows a closure model. However, for the near 
wall region in boundary layers expressions can be given that are exact in the wall limit. 

- Avd 2 {v\v' 2 )dl 2 {v 1 ) (83) 

S 22 = S 3 5 3 = 0.0 (84) 

S 5 12 = - ud 2 (v' 2 2 )d 2 22 ( Vl ) (85) 

The model assumption is essentailly the assumed validity of these expressions for finite dis- 
tance from the wall. Since all expressions are proportional to the laminar viscosity, quick 
decay with increasing wall distance can be expected. 
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2.0 Reynolds stress and complete second order models. 


The usual closure models based on the trace of the dissipation rate tensor will be considered 
first. Several versions of the Reynolds stress model are available and the most relaible stress 
model (Vandromme et al., 1983) will be evaluated. It will serve as test bed for the complete 
second order closure. 


2.1 Stress equations. 

The balance for the Reynolds stress components is given by 

d , d , , , , w , , ,,d(vff) . , , _ Jx d{vg) 


+ V(f£ + g» 

a , $<«>» 


+ 


<9.r 


•(/'- 


<9x. 


- (P)( v 'a v 'l) v 'y) ~ (P'»i) “ tpyip'v'c,)) ~ (/>)« mJ 


( 86 ) 


, «y W 

where the dissipation rate tensor is defined by (1). The split of the pressure correlation into 
pressure transport and pressure rate of strain correlations can be shown to be inappropriate 
near the wall. It follows from (19) that the Taylor series for the pressure gradient is given in 
terms of velocity derivatives at the wall up to second order. Hence, it is possible to represent 
the Taylor series for correlations involving the pressure gradient in terms of local velocity 
correlations at the wall up to second order. If the pressure correlations are split as in (86) 
the pressure fluctuation itself appears and the solution of the Poisson equation introduces 
the well known integral contributions. 

The standard second order closures employ the equation for the trace of the dissipation 
rate tensor 

6 — 7}{ e aa) 

which follows at once from (2) 

(d t + {vjd y )e = d y [ud y e - (v\e) + F?} + S 1 + S 2 + S 3 + S 5 - D (87) 

The source terms are defined as follows. There are two groups of secondary production terms 

S 1 = -{e ay )d 1 {v Q ) (88) 

and 


S 2 = -(vdf,v' Q d y v' a t )(d y (vs) + d6{v y )) 
The primary production term is given by 

S = —{^a-,dyV Q ) 


(89) 


(90) 
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and the pressure correlations can be given in the form 


5 4 = --(^p'^0 (91) 

P 

It was shown in chapter 1.2.1 that the pressure correlations can be represented as the sum of 
the divergence of a flux and a redistributive source which has zero trace. Hence we get 


S 4 = Os F? 


(92) 


where the flux is defined by 


or 


2v 

Ffr — ” £u6i)€ury<> dijV Q ) 


Fl = -j(a,p'a,v' s ) 


(93) 


(94) 


The near wall production of dissipation rate is proportional to the curvature of the mean 
velocity profile 

S 5 = -2v(v' 6 dy Q )0lt(v a ) (95) 

and the viscous destruction of dissipation rate is given by 


D EE 2^X$X> 


(96) 


The properties of the source terms have been established in chapter 1. and we observe that 
none of the source terms is closed in contrast to the equations for the dissipation rate tensor. 


2.2 Standard second order closure model. 


The development of the full second order closure was based in a systematic way on 
existing closure schemes. The standard second order closure using the trace of the dissipation 
rate tensor set up by Vandromme et al. (1983) was an important stepping stone and can be 
applied to test closure models for the various processes governing the dissipation rate tensor. 
The model is given by 


W<! + < ><**> = + Wifr) 




Ox. 


Ox. 


• ^ C/ ll 

The pressure rate of strain correlations are modelled by 

- €,/„(/?* )(/>)|((y>^) - p afi k) - f P (R e )(p){^y^(P Ql} - 1 6 a0 P) 


( 97 ) 
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( 98 ) 


Scz-2 2 30c, - 2 d{v s ) 9{v„) 

— — {Dag ~ 3 KsP) ~ 55 + l^T 1 + **' } 

as the sum of return to isotropy, fast response and wall contributions, where 


and P = l/ 2 P ao and 


r - - - 1«?7 

■ -<-X>« - wig 


(99) 


( 100 ) 


The wall contribution is given by 

e 

k 


K>=M{c\ t.(WM ~ f «<**) + 4(^.0 - A*) + 

3 dx y ' ej/ 


(101) 


The dissipation rate tensor must be modelled in terms of its trace and the anisotropy of the 
stress tensor 

+ (i - /.)§«.**) (102) 

where f 3 depends on the local Reynolds number. The turbulent flux is modelled by 


k. , d 


Fal3-r= ~ C»t(p)-{v y V s ) — (v a V 0 ) 


(104) 


The constants are given by c 3t = 0.25, c x = 1.5, c 2 = 0.4, c\ = 0.1597, c' 2 = 0.0133, 
C 5 = 0.0041. The low Reynolds number functions are given by 


fs(Re) = exp 


(i + P e /3o) 


and 


f p (R e ) = tanh(R e / 50) 


(104) 


(105) 


The equation for the trace of the dissiaption rate tensor requires closure for all relevant 
processes on the right hand side. Following Vandromme et al. (1983) the model 

(# + + ‘■‘wl-l - 


dx y 1 dx 


dx. 


« ( djk\\ n r . , , . € d 2 {v Q ) 3 2 (v a ) 
-c, 2 (p)MR,)jl' - 2*(— j ] + 2 


dx a 
(106) 
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The low Reynolds number functions are given by 


/*' = [!- exp(-/ l/1 i? s )] 2 (l + &) 


(107) 


and 


/ i = l 


1 

8 (/„ + 10 -»®)* 


(106) 


where 


/ 2 = l-0.22exp(-S) 

0 


Ry 


V 


(109) 


(HO) 


and = 0.0165, f „ 2 = 20.5. The performance of this closure model was tested in flat plate 
boundary layers. The results are contained in fig.l to fig.8. The mean velocity in fig.l is in 
good agreement with the law of the wall. The turbulent Reynolds number and the damping 
function f,(R e ) in fig.2 and fig.3 prove that the boundary layer is fully developed. The 
Reynolds stress components in fig. 4 to fig.8 show the expected distributions with maximal 
values in reasonable agreement with the experiments. 



Figure 1 Variation of normalized mean velocity 
with wall distance in wall units 
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Six Equation 
Model 

-2 : x/D Location 2.013 

: x/D Location 4.328 

-6- : x/D Location 6.655 

—8“ : x/D Location 9.014 



Figure 2 Variation of turbulent Reynolds number(R ex) 
vs. normal distance from wall in wall units 



Figure 3 Damping function (F s )for dissipation rate 
term in Reynolds stress equation vs. y + 
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Six Equation 
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: x/D Location 2.013 
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-■8- 
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Figure 4 Variation of normalized fluctuating velocity u 
with wall distance in wall units 
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Figure 5 Variation of normalized fluctuation velocity 
with wall distance in wall units 







Figure 6 Variation of normalized fluctuating velocity w' 
with wall distance in wall units 



Figure 7« Variation of normalized sheer stress u'v' with 
wall distance in boundary layer thickness 
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i in • f f 

Figure 8 Variation of normalized turbulent kinetic energy 
with wall distance in wall units. 
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2.3 Full second order closure model. 


The development of the full second order closure was prepared in the first chapter which 
contains the properties of all transport equations for second order moments derived from 
series expansions with respect to wall distance. The results obtained with series expansions 
are called growth laws since they describe in first or higher order the growth of correlations 
with wall distance. The challenge is now to construct a closure model that satisfies all tensorial 
and realizability conditions and the growth laws. 


2.3.1 Closure model for the stress equations. 


The closure model for the stress equations follows closely the model developed by Van- 
dromme et al. (19S3). The only difference is that the model (102) for the components of the 
dissipation rate tensor is not used. The stress equations appear then according to chapter 
2.2 as 


W(f + w£)(»M> = 


-ci/ P (i2e)(p)|((t«)-!^)-/ p (/2 e )(p){ n 


c 2 + 8/n 2 r 802-2,^ 2 


■(P a! 3--S a pP)- 11 


( D a0 - -SaijP) 


30cj_— _2 d(v/j) 
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dx. 




- (p)M 


(111) 


where P Q p and D Qi j are defined in (99) and (100) respectively and the wall contribution 
to the pressure correlations is given by (101). The low Reynolds number functions are all 
established in chapter 2.2. 


2.3.2 Closure model for the dissipation rate equations. 


The transport equation for the dissipation rate tensor 


d , . d . , , dF Q p 7 . - , . , 

-^ + Sit + Sig + Sh, + si, + Sif - Da, (112) 


l dt )('■.<») - 


requires closure for the turbulent flux contained in the total flux F„ ^ and the source terms 
$00' $0)3 an< ^ Dafl- 
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Model for the turbulent flux: This model follows the suggestion of Vandromme et al. 
(1983) for the turbulent flux of the trace of the dissipation rate tensor 


■F 1 a()~i — v 


d[tqp) 

dx 1 


+ c f{p) 


k 





djeap) 

dx& 


(113) 


The only modification is the use of the time scale (69) to improve the behaviour neax the 
wall. The value for the constant c ( =0.1 is consistent with the six-equation model. A more 
sophisticated model was developed in chapter 1.3.2. 

Model for the source 5^: The secondary production can be neglected near the wall 
according to chapter 1.2. 2. 4. 

Model for the source S The primary production is not of leading order near the wall 
but varies like the closed production term S l afi in the near wall region according to ch.1.2.2.5. 
The present model utilizes this property and the analogy to the stress transport exploited 
in the equation for the trace (106), where the high Re-number part of the model for the 
difference between the leading terms having the character of a production term is modelled 
proportional to the production of kinetic ernergy 


o3 ~ _ /.w/_ d ( v °) , - 

s °>- ~ c<3 + '“•'irr 1 


c.jaR',)(p) f ( n 4 ) 


dx^ 


dx. 


The constants are given by c (1 = 1.45 and c {3 = 1.0. 

Model for the source term The pressure correlations were shown to be of leading 

order near the wall for all dissipation rate components except €22 (see ch.1.2.2.9). They can 
be split into diffusive and reditributive source terms (ch. 1.2.1). The diffusive part is assumed 
to be represented by the closure for the turbulent flux (113). The model for the redistributive 
source is analoguous to the return to isotropy model for the stress transport equations. It is 
given by 


S Q p— C f 4 {p) 


(115) 


with c e4 = 12.5. 

Model for the destructive term D a p : It was shown in chapter 1.3.4 that the high Re- 
number model for the difference of primary production and viscous destruction becomes 
singular as the wall is approached. Furthermore, the growth law for the viscous destruction 
term D a p implies that a model similar to the high Re-number case 


D a p=c (2 (p) 


(fop) 

T 


with a scalar time scale is impossible. It follows that a tensorial time scale m ust be con- 
structed to conform with the growth law. The present closure is a composite expression 
containing the high and the low Re-number models 


D Q p=c ( 2 \{p)(\ — fs(Ret)) 


( k ) 
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djtgj) djepj) 

dx v dx $ 
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+c ( 22 (p)n Q n li n y n 6 ( 1 - f a (R et ))j(e^&) + c (23 (p)/ a (i? et ) n ' r?? ^ €7 ^ (e a0 ) (116) 

where the preliminary values for the constants are c ( 21 = 0.32 • 10~ 4 , c c 22 = 26.25, and 
c f 23 = 12.5. The low Re-number functions are set up as follows: 

90.5 

fu{Rtt) = (1 - exp[— 0.0165 -Rj,]) 2 (1 - z =— ) (117) 

K e t 

where 


denotes the dimensionless wall distance and 

fs(Ret) = tanh(0M4R £t ) (118) 



The turbulent Reynolds number is defined by 


R 


et. 


£ 

ev 


2.3.3 Preliminary results for the complete second order closure. 

The system of nine parabolic differential equations was tested in reduced form by pre- 
scribing the profiles for mean velocity and the Reynolds stress components which were ob- 
tained with the six equation model discussed in chapter 2.2. The numerical solution for the 
remaining equations for the disssipation rate tensor was carried out and convergence was 
achieved after a few hundred steps. The results are presented in fig. 9 to fig. 16. The figures 
contain also as broken line the dissipation rate components deduced from the local relation 
suggested by Launder and Reynolds (1983) and modified by Lai and So (1990) 


{^ap ) — jj (1 fw(Ret))^Qf}^ 


, , , D , e {v'aVfi) + {v’av'Jn^n,, + (vj,i >' 7 )n 7 n tt + n Q n/?n 7 n«(v;^) 

+ /u)(-Uet)— ~ - — ~ ~ — (119) 


where 


1 + Jk n y n s( V 'A) 


/„(«„) = exp[-(^) 2 | 


The dissipation rates in the figures are normalized with the wall variables v and u T 

(w)" 


e + — 
€ a/3 ~ 


U, 
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with u T = \Jt w / p denoting the wall shear velocity. The symbols in the figures represent the 
results of the numerical experiments carried out by Mansour, Kim and Moin (1988). 

The component in fig. 9 (logarithmic scale) and fig. 10 (linear scale) shows a negative 

gradient near the wall which implies that the time scale 


T 


-1 


1 0 - ( MV 

k k\dy ) 


would change sign and the destructive term would become a production term violating real- 
izability. The prediction of in the outer (fully turbulent) part of the flow field is too small 
and the near wall part appears too large compared to the direct simulation, but the profile 
shape is in good agreement with the numerical experiment. The component e 22 in fig. 11 and 
fig. 12 reasonable agreement between the full second order closure and the direct simulation, 
but the local relation of Launder and Reynolds overpredicts the component by a factor of 
three. The prediction of the component in fig. 13 and fig.14 shows a similar behaviour as 
the component ejj in fig. 9/10. The shear dissipation e+ 2 in fig. 15/16 shows overprediction by 
the full second order closure and underprediction by the local relation near the wall whereas 
the outer part is in good agreement with the direct simulation data. 


2.3.4 Conclusions. 

The results presented lead to several conclusions. It is clear from the theoretical develop- 
ment that only the full second order closure offers the tools to construct the appropriate time 
scales in the near wall region of a turbulent boundary layer. The growth laws for the various 
correlations appearing in the stress and dissipation rate balances limit severely the model 
expressions and indicate that composite models for the high and the low Reynolds number 
regimes must, be established. The model discussed in this chapter produces good results if 
velocity and stress components are held fixed and this indicates that the model expressions 
are consistent with the direct simulations. However, the stability of the closure needs to be 
investigated and this part may lead to modifications of the present version of the full second 
order closure. This part of the project is currently under way. 
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Figure 9 Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Mansour, 
Kim, and Moin direct simulation results. A Launder & Reynolds. 
— : Present Work. 



Figure 10 : Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Mansour, 
Kim, and Moin direct simulation results. A Launder & Reynolds. 
— : Present Work. 

s 
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Figure n : Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Mansour, 
Kim, and Moin direct simulation results. A Launder & Reynolds. 
— : Present Work. 



Figure 12 : Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Mansour, 
Kim, and Moin direct simulation results. ▲ Launder & Reynolds. 
— : Present Work. 
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Figure 13 : Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Man sour, 
Kim, and Moin direct simulation results. A Launder & Reynolds. 
— : Present Work. 
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Figure 14 : Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Mansour, 
Kim, and Moin direct simulation results. ▲ Launder & Reynolds. 
— : Present Work. 
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:Prcsent Work. 



Figure 16 .'Variation of Psuedo 9 equation model of normalized 
dissipation rated with wall distance. O : Results from Mansour, 
Kim, and Moin direct simulation results. A Launder & Reynolds. 
— : Present Work. 
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